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Shear-strain and shear-stress correlations in isotropic elastic bodies are investigated both theoreti¬ 
cally and numerically at either imposed mean shear-stress r (A = 0) or shear-strain 7 (A = 1) and for 
more general values of a dimensionless parameter A characterizing the generalized Gaussian ensem¬ 
ble. It allows to tune the strain fluctuations /i 77 = fiV (Sj 2 ) = (1 — A)/G eq with /3 being the inverse 
temperature, V the volume, 7 the instantaneous strain and G eq the equilibrium shear modulus. Fo¬ 
cusing on spring networks in two dimensions we show, e.g., for the stress fluctuations fi TT = fdV(5f 2 ) 

(f being the instantaneous stress) that /i TT | A = /a\ — AG eq with /jla = y TT \\=o being the affine shear- 
elasticity. For the stress autocorrelation function C TT {t) = f3V (5f(t)5r( 0)} this result is then seen 
(assuming a sufficiently slow shear-stress barostat) to generalize to C TT (t )| A = G{t) — AG eq with 
G(t) = C TT (t) |a=o being the shear-stress relaxation modulus. 

PACS numbers: 05.70.-a,05.20.Gg,05.10.Ln,65.20.-w 


INTRODUCTION 


Background. A “simple average” A = (A) of an ob¬ 
servable A pQ does not depend on which thermodynamic 
ensemble it is measured in, at least not if the system 
is sufficiently large and if each ensemble samples indeed 
the same thermodynamic state point DM- An exam¬ 
ple for such a simple average is the affine shear-elasticity 
/i a characterizing the energy change of an affinely shear- 
strained elastic body I3HD1 as properly defined below in 
sect. Ill As one may verify numerically [ 8 ], it is irrele¬ 
vant whether one computes fi a in the NVyT-ensemble 
at constant particle number TV, volume V, shear strain 7 
and temperature T = 1//3 or in the conjugated NVtT- 
ensemble where the strain is allowed to fluctuate sub¬ 
jected to the constraint that the internal mean shear 
stress r is imposed by the external applied shear stress 
^ext 0 . In contrast to this, the fluctuation ( 5A5B) of 
two observables A and B may depend on whether an ex¬ 
tensive variable X or its conjugated intensive variable I 
is imposed DM- Focusing in the present work on shear- 
strained isotropic elastic networks, as sketched in panel 
(a) of fig. [l] the relevant extensive variable is the rescaled 
shear strain X = V 7 and the conjugated intensive vari¬ 
able the shear stress I — r. Using the Lebowitz-Percus- 
Verlet (LPV) transform [ 2 ] it is seen for the (rescaled) 
mean-squared fluctuation (i TT = f3V (Sr 2 ) of the instan¬ 
taneous shear stress f that 00 


Rtt |a=0 — Rtt |a=1 “I - 


eq 


(1) 


with G eq being the equilibrium shear modulus. For later 
convenience the NYrT-ensemble is indicated by “A = 0” 
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FIG. 1: Sketch of addressed problem: (a) Shear strain ex¬ 
periment in the {x, y) -plane with shear strain 7 and applied 
external shear stress r ex t- (b) External potential I/ ex t{p()/V 
according to eq. § with the bold line indicating the external 
spring of spring constant G ex t- (c) Dimensionless parameter 
A = G ex t/(G eq + G ex t) comparing the external spring con¬ 
stant G ex t and the equilibrium shear modulus G eq - We set 
G eq = 16.3 as for the elastic model system considered below. 


and the NVyT-ensemble by “A = 1”. Generalizing eq. 0 
in the time-domain it has been shown for the stress-stress 
correlation function C rr (t) = /3V(5r(t)5f( 0 )) that [9lflQ] 

C T r{t) |a=o — C rr (t) |a=i + G eq . ( 2 ) 

This assumes that the shear-barostat needed to sample 
the NVrT-ensemble must act very slowly on the time- 
scales used to probe C TT (t) |a=o- This may be realized 
equivalently by averaging over an ensemble of configu¬ 
rations with quenched strains 7 distributed according to 
the NVrT-ensemble 0. It is due to the ergodicity break¬ 
ing generated by the averaging over the quenched ensem- 
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ble that C rr (t) |a=o G eq stays finite for t oo, while 
C TT (t) |a=i must vanish for large times (as one commonly 
expects for all correlation functions in ergodic systems 
[IT]). By integration by parts it is seen that C TT (t) |a=o 
is equivalent to the experimentally important shear re¬ 
laxation modulus G(t) The transform eq. © thus 
implies G{t) = G eq + C rr (t) |a=i, he. the relaxation mod¬ 
ulus may be measured using the equilibrium modulus G eq 
and the correlation function C tt (£)|a=i both determined 
in the NVyT-ensemble [9j. 

Generalized Gaussian ensemble. As sketched in panel 
(b) of fig. [l] it is straightforward to interpolate between 
the NVyT-ensemble and the NVrT-ensemble by impos¬ 
ing an external field 

G e xt (q)/h^ — ^”ext('T Text) T~ Text) (3) 

with G ext being the spring constant of the external har¬ 
monic spring. The standard NVrT-ensemble at r = r ext 
is recovered by setting G ext = 0 . Note that our approach 
is conceptually similar to the so-called “Gaussian ensem¬ 
ble” proposed some years ago by Hetherington and others 
nans] generalizing the Boltzmann weight of the canon¬ 
ical ensemble by an exponential factor U ext (E) oc E 2 of 
the instantaneous energy E. A similar external spring po¬ 
tential has been also used in the recent “elastic bath” ap¬ 
proach by Workum and de Pablo [14j . Choosing the refer¬ 
ence strain 7 ext equal to the mean strain 70 of the NVrT- 
system at vanishing shear stress r = (f) = r ext = 0 allows 
to work at constant zero mean shear stress irrespective 
of the strength of the external potential [IT. All the 
ensembles considered correspond thus to the same ther¬ 
modynamic state, i.e. all first derivatives of the energy or 
the free energy [3] and all simple averages are identical. 
As sketched in panel (c) of fig. [l] G ex t is not necessarily 
positive, reducing the strain fluctuations, but may even 
become negative, which thus amplifies the fluctuations. 
It has been argued [14j that this may allow a more con¬ 
venient determination of the elastic modulus. Since the 
external spring is parallel to the system, the combined 
system and external device have an effective spring con¬ 
stant G e ff = G eq + G ext . Defining the dimensionless pa¬ 
rameter A = G ex t/(G eq + Gext) the strain fluctuations of 
the combined system are thus given by 


M 77 = pv (5f) = 1/Geff = (1 - A)/G eq , (4) 

i.e. must vanish linearly with A. NVrT-ensemble statis¬ 
tics is expected for A 0, while NVyT-statistics should 
become relevant in the opposite limit for G ex t oo, he. 
A 1. The system must become unstable in the limit 
G ex t ~^ — G eq , i.e. A —>• —oo. 

Key results. The aim of the present study is to gen¬ 
eralize the relations for static fluctuations, eq. ([!]), and 
dynamical correlation functions, eq. (J2|, to the more gen¬ 
eral transformations between Gaussian ensembles char¬ 
acterized by the continuous parameter A < 1. In this 
way we want to make manifest that these transformation 


relations are generated by the continuous change of the 
constraint imposed on the fluctuations of the extensive 
variable. Focusing on spring networks in two dimensions 
we show, e.g., for the stress-stress fluctuations fi TT in dif¬ 
ferent A ensembles that 

hrr | A = hrr |a= 0 ^G eq (5) 

with fi TT |a=o being given by the affine shear-elasticity /i a 
mentioned above. Assuming a very slowly acting shear- 
barostat, which is irrelevant for the system evolution for 
short times t r*(A), the above result is then seen to 
generalize in the time-domain for the stress-stress corre¬ 
lation function 


C TT (t) |a = C TT (t) |a=o - AG eq for t < T*( A). (6) 

Since G(t) = C TT (t) |a=o, eq. © allows the determina¬ 
tion of the relaxation modulus G(t) from the equilibrium 
modulus G eq and the stress-stress correlation function 
C TT (t )|a for any A and t <C t*(A). The upper time limit 
t*(A) is seen to be set by the diffusion time T 77 of the 
instantaneous strain j(t) over the typical strain fluctua¬ 
tions for A. As a consequence from eq. © and eq. © 
it appears that it is the equilibrium modulus G eq of the 
system, which generates both transforms 

Jj MttL = A C TT {t)\ x = -G eq (7) 


with fjb TT |a = /a a and C TT (t) | A = G(t) for A = 0. Similar 
relations, allowing also to determine the creep compliance 
J(t) [HHZ], are obtained for strain-strain and strain- 
stress correlations functions. 

Outline. The stated key relations eq. © and eq. © 
are justified theoretically in sect. [TTJ The static fluctua¬ 
tions are discussed in sect. |II A before we address the dy¬ 
namical correlation functions in sects. pTBlITT Clan d |!TDl 
and the macroscopic linear response in sect. 11 Algo¬ 
rithmic details are given in sect. Ill where the specific 
model system considered is introduced in sect. Ill A[ The 
canonical affine plane shear transformations used in this 
work are specified in sect. IIIB[ the instantaneous shear 
stress f and the instantaneous affine shear-elasticity jl\ 
in sect. |III C| The zero-temperature ground state prop¬ 
erties of the two-dimensional elastic network are summa¬ 
rized in sect. nma Some technical details related to the 
finite-temperature simulation of the Gaussian A-ensemble 
using a Metropolis Monte Carlo (MC) scheme as a func¬ 
tion of the maximum strain displacement rate ft, the sec¬ 
ond key operational parameter of this study, are given in 
sect. ITTTeI Section [TV| presents the numerical results ob¬ 
tained by molecular dynamics (MD) simulations at one 
finite (albeit small) temperature T and one relatively 
high value of the friction constant £ of the Langevin 
thermostat used [l]. The relevant static properties are 
described in sect. |IV A before we turn to the dynam ical 
strain-strain (sect . |IVB ), strain-stress (sect. |IV C ) and 
stress-stress (sect. |iVDf correlation functions. Our work 
is summarized in sect. m 

































3 


II. THEORETICAL CONSIDERATIONS 
A. Static fluctuations 


Lebowitz-Percus-Verlet transform. We begin by re¬ 
stating the LPV transform in a convenient form assuming 
that the relevant extensive variable is the (rescaled) shear 
strain X = V7 and the conjugated intensive variable the 
shear stress I = r. Following Lebowitz et al. mu , one 
verifies (see also sect. II.A of ref. [8]) that to leading 
order 


6 MB) 


= (§Asb) 

A=0 \ / 


G e q c)A c)B 

X=1 flV dr dr { } 


for the transformation of ( SASB) with 5 A = A — (A) 
and SB = B — (B) from the NVyT-ensemble (A = 1 ) to 
the NVrT-ensemble (A = 0 ). The more general transfor¬ 
mation between arbitrary A-ensembles can be found by 
reworking the saddle-point approximation of ref. [2] tak¬ 
ing into account that the fluctuations around the peak of 
the distribution of the extensive variable is now not set 
by the modulus G eq of the system but by the effective 
modulus G e ff = G eq + G ex t- How this may be done can 
be seen in sect. 2.5 of ref. m for the volume X = V as 
extensive variable and the (negative) pressure I = —P 
as intensive variable. Rewriting the latter result to the 
present case we get the generalized LPV transform 


sAsb) 


x = (sAsb) 


+ ( 1 ~ x )wPir (9) 

A=1 pV OT OT 


which reduces for A = 0 to eq. ©• 

Strain-strain fluctuations. Let us check the ensem¬ 
ble dependence of the rescaled strain-strain fluctuations 
/i 77 = flV (£y 2 ) p~ 9 ] . The generalized Gaussian ensem¬ 
ble corresponds to replacing the shear modulus G eq of 
the system by the effective modulus G e ff = G eq + G ex t- 
As already stated above, eq. © , this leads to pQ] 


/i 77 — 1 /G e ff — (1 — A)/G eq . (10) 


Stress-stress fluctuations. We turn to the most im¬ 
portant stress-stress fluctuation fi TT = flV (dr 2 ) [T 9 ] . We 
set now A = B = r in the LPV transform. Since the 
stress fluctuations do not vanish in the NVyT-ensemble 
the corresponding term must now be included. This 
yields [20] 


t^TT |a — !^tt\\=i T (1 A)G eq . ( 12 ) 


Interestingly, the contribution (1 — A)G eq may be rewrit¬ 
ten using the notation (/(q)) 7 = f dy p(y; A)/(7) for 
the strain-average of a property /(y) with p(y; A) be¬ 
ing the normalized strain-distribution for the consid¬ 
ered A-ensemble. The total mean stress r of the en¬ 
semble is thus given by r = (r(y)) 7 with r(y) being 
the average shear-stress of all configurations of shear- 
strain 7. Using that p(y; A) is a Gaussian and that 
Sr{fl) = r(y) — r = G eq (y — 70) with 70 being the maxi¬ 
mum of the distribution, it is then seen that 

(l-A)G eq =|| = /3Fqr(7) 2 ) 7 . (13) 

Using eq. ( p~2| ) this implies in turn 

MttIa = Mtt|a=i + PV (St( 7) 2 ) t (14) 


as one expects to lowest order from a standard saddle- 
point approximation. 

Stress-fluctuation formula for shear modulus. Sub- 
stracting the transform for A = 0 from eq. (12) confirms 
the key result eq. (§. As shown in ref. m , h'rr |a=o can 
be reduced by integration by parts to the affine shear- 
elasticity /iA- Since the latter expression is a simple av¬ 
erage, i.e. the same value is obtained in any ensemble, 
eq. © may be further simplified as m 

hrr |a h'A AG eq . ( 15 ) 


It is thus sufficient to compute the material constants 
/a a and G eq in any ensemble to obtain the stress-stress 
fluctuations fi TT as a function of A. Note that the spe¬ 
cial case for A = 1 corresponds to the well-known stress- 
fluctuation formula gnu nu [22] 

G e q = /a a AUt|a=i (15) 


One verifies readily that the postulated LPV transform 
for general A, eq. ([ 9 |, is consistent with this result. To 
see this one sets A = B = V7 and uses the fact that the 
strain fluctuations must vanish in the NVyT-ensemble, 
i.e. that the first term on the right hand-side of eq. © 
must vanish. 

Strain-stress fluctuations. Setting A = V7 and B = r 
and using again that the NVyT-term in eq. © must 
vanish, it is seen from the LPV transform that the strain- 
stress fluctuations fi 1T = flV (d^/Sr) [ 19 ] should scale as 

fJbjr = 1 — A. ( 11 ) 


This result can be also obtained directly by replacing in 
the definition of n 1T the fluctuation Sf by G eq ^7 and 
using then the strain-strain relation, eq. (10). 


used in numerous numerical studies to compute the 
modulus G eq conveniently in the NV7T-ensemble 0- 

mmmm- 

B. Dynamics: Definitions and general relations 

Introduction. We define now several dynamical ob¬ 
servables of interest, which will be investigated numer¬ 
ically in sect. |IV[ and remind some well-known gen¬ 
eral relations DmnEzi. Time translational symmetry 
+ St) and time reversal symmetry (t —t) are as¬ 
sumed. We use a(t) and b(t) for either the instantaneous 
shear strain 7 or the shear stress f, a = (a) and b = (b) 
for their thermodynamic averages and Sa(t) = a(t) — a 
and Sb(t) = b{t) — b for their time-dependent fluctuations. 
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Mean-square displacements. We define (generalized) 
mean-square displacements (MSD) by 

5ab(i) = ((“(*) - “(°))(K*) - K°))) ( 17 ) 

where the prefactor /3V /2 has been introduced for conve- 
nience [ 32 ]. Obviously, g ah (t) = gb a (t), g ah (t) = g a b(~t) 
as may be seen using the above-mentioned symmetries, 
9ab if) 0 for t 0 and 


9ab(t) Mab = PV ( 6a 6b ) for t > T ab (18) 

with T a b being the characteristic time needed to reach 
this thermodynamic limit [28] . 

Correlation functions. Similarly, we define dynamic 
correlation functions by C a b(t) = /3V (Sa(t)Sb( 0 )) [19]. 
Obviously, C ab (0) = g ah , C ab (f) = C ab (-f) and C ab (f) = 
C ba (f). The latter identity is seen from [3l HT] 

(a(t)b( 0)) = (a(0)6(-i)) = (&(f)a(0)) (19) 

where the time translational invariance is used in the first 
step and the time reversal symmetry in the second step. 
Using again both symmetries one verifies that nn Ha 

9ab if) = Cab(0) - C a b(t) = 9 ab C a b(t). (20) 

This implies that C a b(t) —>► 0 for large times t T ab [28] . 
Albeit g a b(t) and C a b(t) thus contain the same informa¬ 
tion it will be sometimes better for theoretical or numer¬ 
ical reasons to focus on either g a b(t) or C a b(£). Since 
g a b(t) 0 for t —)> 0 it will be natural, e.g., to consider 
the double-logarithmic representation of g a b{t) to clarify 
the power-law scaling of the correlations at short times. 

Finite-time dependence of fluctuation estimation. As 
defined in sect. EE fi a b is a thermodynamic ensemble 
average, hence, a time-independent property. In practice, 
fi a b may, however, often be estimated using a time series 
(afe, b k ) with a finite number n of more or less correlated 
entries mis]. It is supposed here that these series are 
sampled with a constant time interval of length St , i.e. 
n corresponds to a time t = (n — 1 )St. With A k m 
n Sfe=i A k denoting such a finite average, one samples 
the n-dependent observable m 


Mab( n ) = PV ( («fc - a k )(h - b k ) 


( 21 ) 


which can be verified by straightforward algebra. (See 
sect. 2.4 of ref. [37].) Using time translational invariance 
this implies 


2 x A 

Kb( n ) = g&b ( 6t ( 23 ) 

k =1 


where the weight (n — k) stems from the finite length of 
the trajectory used. We change now to continuous time 
variables, k —> s = (k — l)5t, and replace the discrete 
sum by the time integral 


2 

Mab(b = ~ t J o ds (1 - s/t) g ah {s). 


Using eq. (|20|) this yields the general relation 


1 - 


Mabft) 

Rab 


_ 2 [* 
~tj 0 


ds (1 — s/t ) 


c ab (s) 

c ab (o) 


(24) 


(25) 


which can be used to obtain /x ab (t) if the correlation func¬ 
tion C a b{t) is known. 

Relaxation time 0 a b. Interestingly, the time integral 
in eq. ( [25] ) must become constant in all cases considered 
below where C a b(t) vanishes ultimately for t oo. (This 
applies if a finite shear-barostat is switched on.) Defining 
the characteristic relaxation time 


0 ab = lim [ ds (1 - s/t) tEFH ; 

Jo Mab 

this leads to 


(26) 


1 - MabW/Mab 2 0 ah /t for t oc (27) 


as one expects for the Poisson statistics of uncorrelated 
events [5]. One may thus determine the relaxation time 
0 a b from the large-time asymptotics of p ab (£) [29] . 


C. Dynamics: Additional model assumptions 

Let us suppose that the MSD g a b{t) is diffusive for 
short times, i.e. g a b(£) = D a bt/2 for t <C T a b with 
D a b being a diffusion constant. Matching this short time 
regime with g a b{t) = fi a b for t T a b gives the possibility 
to operationally define T a b as the crossover time 


with (...) standing for an additional ensemble average 
over different time series of n subsequent data points. We 
have / 2 ab (n) = 0 for n = 1 and p ab (n) —>* Mab for n —)> oo 
in an ergodic system. Interestingly, the detailed time- 
dependence of b (n) can be obtained from a weighted 
integral of the correlation function C a b(t) |5j E- To see 
this we note first the identity 

-;-— i n 

(dk - a k )(b k - b k ) = ^2 ^2 ~ ~ bi) (22) 

k,l =1 


2~ab — 2/i a b/ D a b. 


(28) 


Different short-time dynamics may suggest of course a 
different operational definition. Let us further assume 
an exponentially decaying correlation function C a b(t) = 
/i a bexp(— x) with x = t/T a b. Using eq. (20) this is seen 
to be consistent with a diffusive short-time regime and 
eq. (|28|). It follows by integration from eq. (25) that 


1 - 


Mabft) 

Mab 


= /Debye(^) = ~^[exp(-x) ~ 1 + x] (29) 
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using the Debye function well-known in polymer science 
m- For large reduced times x 1 this leads to 1 — 
~p a b(t) /Mab 2T a b /£, i.e. by comparison with eq. (27) we 
have T a b = #ab for an exponentially decaying correlation 
function. 

As we shall see in sect. IV D[ it may also occur that 
the correlation function is more or less constant between 
a local (barostat independent time) ta up to a (barostat 
dependent) time r*, i.e. essentially C a b(t) ~ ciL(r* — t ) 
with c being a constant and H ( x ) the Heaviside function. 
It follows from eq. (26) that 


Mab(b _ f C for T A « t « r* 

Mab ~\cr*/£ for 


(30) 


i.e. # a b ~ r* c/2 to leading order. 


D. Dynamics: Ensemble effects 

Introduction. We have just stated various general re¬ 
lations applying to all ensembles. These relations do, 
however, not allow to relate the dynamical correlations 
between different A. Since some ensembles are more read¬ 
ily computed than others, it would be useful to have a 
transformation relation such as the LPV transform for 
static fluctuations considered in sect. EE We emphasize 
that the correlation functions depend in general on the 
ensemble, i.e. the parameter A, and on the dynamics of 
the shear-barostat. Since in the end we want to describe 
the intrinsic relaxation dynamics of the system, it is suf¬ 
ficient to focus on the limit where the barostat becomes 
very slow such that it becomes essentially irrelevant for 
the system evolution below an upper time t*. We con¬ 
sider the limit where this time r* is be much larger than 
any intrinsic relaxation time of the system. Since we want 
still to consider meaningful thermal averages with respect 
to the chosen ensemble, we need either trajectories of a 
time interval t tra j much larger than r* to allow the full 
sampling of the phase space or, equivalently, we need to 
average over independent start configurations represent¬ 
ing the ensemble. Under these two assumptions useful 
transformation relations can be formulated by reworking 
the generalization of the LPV transform for dynamical 
correlations replacing the system modulus G eq by the ef¬ 
fective modulus G e ff. Being beyond the scope of this 
paper, we proceed by postulating the central scaling re¬ 
lation for the MSD g a b{t) and argue briefly why this rela¬ 
tion is natural. We discuss then in turn the strain-strain 
correlations G 77 (£), the strain-stress correlations C 1T (t ) 
and the stress-stress correlations C rr (t) for different A. 
Alternative, more direct ways to derive the relations are 
mentioned. 

Scaling relations. We postulate that under the two 
assumptions made above the MSD g a b(t) does not depend 
on the ensemble 


(31) 


i.e. the MSD behaves as a simple mean. This scaling 
postulate is justified by two facts. Firstly, the barostat is 
(by construction) too weak to change for the evo¬ 

lution of the system. Secondly, that the starting points of 
the trajectory at t = 0 are distributed according to the 
considered ensemble must become an irrelevant higher 
order effect (vanishing rapidly with the system volume 
V), since the MSD probes displacements with respect to 
the starting points, not their absolute values. The fun¬ 
damental scaling, eq. implies using eq. ( [20| ) that 

G a b (t) = Mab(A) - g a b(t) for t < T*( A) (32) 

where only the first term on the right hand-side depends 
on A. This leads us finally to the general transformation 
relation between correlation functions 


Gab (t)\\ = Gab (t) | A=1 + Mab | A ~ M ab | A=1 (33) 


where A = 1 stands for the NVyT-ensemble with 7 = 7 ext 
The two static contributions on the right hand-side can 
be further simplified using results from sect. |II A| This 
demonstrates the linear relation /i a bU — Mab|A=i ^ 1 — A. 

Strain-strain correlations. Since q(t) ~ 7(0) for t <C 
t*, the MSD g 11 (t) must vanish to leading order. Using 
eq. ( [32] ) this implies [20] 

G 77 (£) = /i 77 = (1 - A)/G eq for t < t*( A). (34) 


This result is directly obtained by setting q(t) = 7(0) 
in the definition of G 77 (t). We emphasize that it is also 
consistent with the LPV transform, eq. ([9|, as one verifies 
by setting / = r, X = V7, A = Vy(£) and B = Uy(0) 
and using that the strain fluctuation term for A = 1 must 
vanish. 

Strain-stress correlations. To determine the strain- 
stress correlation function G 7T (£) one may use again that 
the MSD gy T (t) must vanish since y(t) ~ 7(0) for K<r*. 
Using eq. ( [32] ) this implies 20 

C jr (t) = /i 7r = 1 — A for t r*(A). (35) 


This result is also obtained from the LPV transform set¬ 
ting A = Vj(t) and B = r(0) and using again that 
(SASB)\x=i = 0. The postulated scaling eq. (31) has 
thus led again to a reasonable result. 

Stress-stress correlations. Interestingly, as one may 
see from the NVyT-ensemble limit considered in mm, 
the stress-stress MSD g TT (t) is not expected to simply 
vanish for t < 7 as in the two previous cases. (The 
instantaneous stress f fluctuates even at a fixed strain 
7.) However, eq. ( |33| ) still holds leading to m 


C TT (t)\x = C TT (t )\ X =1 + (1 - A)G eq for t < T*( A) (36) 

where the static term on the right hand-side has been 
simplified using eq. ([5|. This confirms the key relation 
eq. © announced in the Introduction. Note that the 
correlation function C TT (t)\\=i in the NVyT-ensemble 
must vanish beyond some local time scale ta which does 


5ab(i) ~ a 0 for t < T*(A) 
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depend on the network considered but, of course, not on 
the shear-barostat. For a sufficiently slow barostat the 
correlation function thus becomes constant 


C T r{t) = (1 - A)G eq 

where we have dropped |a- 
the remarkable relation 


for r a <C t t*( A) 


Using eq. (32) this 


(37) 
leads to 


g TT {t) = Ma — G eq for r A <S t « r*(A) (38) 


which must hold for all A < 1. We note finally that 
while for A < 1 the MSD g TT (t) — )> fi TT for t t*, 
g TT (t) = /iA ~ G eq holds for all times £ ta for A = 1 
according to stress-fluctuation formula, eq. (16). 


E. Dynamics: Macroscopic linear response 

Introduction. The experimentally important macro¬ 
scopic linear response measured by the creep compliance 
J(t) and the shear relaxation modulus G(t) |T61 fT7l 127] 
may be obtained conveniently in an equilibrium simula¬ 
tion at a given A using some of the correlation functions 
discussed above. Please note that being material prop¬ 
erties of the given state point (experimentally obtained 
using a simple average, not a fluctuation) both response 
functions J(t) and G{t) do, of course, not depend on A. 

Creep compliance. Let us first consider the creep com¬ 
pliance J(t) = (£ 7 (£)) /Sr ext for t > 0 . It is assumed that 
for t < 0 the system is at thermal equilibrium and the 
internal mean stress r equals the applied external stress 
Text of the NVrT-ensemble. After imposing at t = 0 a 
small increment Ar ext , the creep compliance J{t) mea¬ 
sures the ensuing average strain increment (5*f(t)). We 
note en passant that the average internal shear stress 
(f(£)) does neither immediately reach the new equilib¬ 
rium value T ext + 5r ext but shows a similar time depen¬ 
dence as the strain. Reworking the arguments put for¬ 
ward by Doi and Edwards, see eq. (3.67) of ref. [27] . it is 
seen that 


Shear relaxation modulus. The shear relaxation mod¬ 
ulus G(t) = (Sf(t)) /#7 may be obtained from the stress 
increment (Sr{t)) for t > 0 after a small step strain with 
| £ 7 1 < 1 has been imposed at time t = 0. It is well 
known that the components of the Fourier transformed 
relaxation modulus G(£), the storage modulus G'(uj) and 
the loss modulus G"(w), are directly measurable in an os¬ 
cillatory shear strain experiment pum. As seen, e.g., 
by eq. (32) in ref. [TO] . it can be demonstrated by inte¬ 
gration by parts that 

G(t) = C TT (t) |a=o for t < t*(A = 0 ), (41) 


i.e. the barostat should be irrelevant on the time scales 
considered. Using the transformation relation between 


different ensembles, eq. (36]), the relaxation modulus may 
equivalently be obtainec. 


for other A according to 


G(t) = G rr (£)| A + AGeq for t < t*(A). (42) 

Since for large times G(t) —» G eq , this implies C TT (t)\\ 

(1 —A)G eq consistently with eq. (37). For the specific case 
A = 1 eq. (|42|) yields 


G(t) — G tt (£)|a=i + G e 


(43) 


which holds for all times t since the barostat becomes 
irrelevant for A —>> 1 . Note that eq. (43) may also be de¬ 


rived directly (without using the transformation relation 
between different ensembles) using Boltzmann’s super¬ 
position principle for an arbitrary strain history and the 
standard fluctuation-dissipation theorem for the after¬ 
effect function HD as shown elsewhere mm- Two im¬ 
mediate consequences of eq. ( |43| ) are (i) that G(t) only 
becomes equivalent to C rr (t)\\= 1 for t > 0 in the liq¬ 
uid limit where (trivially) G eq = 0 and (ii) that the 
shear modulus G eq is only probed by G(t) on time scales 
t ta where C TT (t) | a=i must vanish. In principle, it is 
thus impossible to obtain the static shear modulus G eq of 
an elastic body only from C TT (t) |a=i as often incorrectly 
assumed mm 


III. ALGORITHMIC DETAILS 


Jit) = fl 77 (i)U =0 for I^Textl <S 1, 


(39) 


A. Model Hamiltonian 


i.e. the creep compliance is most readily computed using 
the strain-strain MSD ^ 77 (t) in the NAGT-ensemble. As 
described in sect. |III E| we shall change the shear-strain 
q(t) using a MC shear-barostat which corresponds to a 
perfectly viscous dynamics. One thus expects G 77 (£) = 
// 77 exp(—£/T 77 ) for the strain-strain correlations. Using 
// 77 = 1 /G eq for A = 0 and eq. ( 20 ) this suggest 


J{t) = — [1 - exp(—£/T 77 )] 


(40) 


To illustrate our key relations we present in sect. IV 


nu¬ 


merical data obtained using a periodic two-dimensional 
(d = 2) network of Ni harmonic springs connecting N 
vertices. The model Hamiltonian is given by the sum 
% = 77id + ^ex of a kinetic energy contribution 


N 

~ m ^ 2 

^id = x y v.j, 


(44) 


with v i being the velocity of vertex i and m its (assumed) 
monodisperse mass, and an excess potential 

Nl 1 

H ex = y] Ui(n) with Ui{r) = -K t (r - Ri) 2 (45) 
1=1 “ 


in agreement with the Kelvin-Voigt model m represent¬ 
ing a purely viscous damper and a purely elastic spring 
connected in parallel. 
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where Ki denotes the spring constant, Ri the reference 
length and 77 = |zj ~ Lj\ the length of spring l. The sum 
runs over all springs l connecting pairs of beads i and j 
with i < j at positions r { and rj . The vertex mass m and 
Boltzmann’s constant &b are set to unity and Lennard- 
Jones (LJ) units p] are assumed. 


B. Canonical affine shear transformations 


While the box volume V is kept constant throughout 
this work, we shall frequently change the shape of the 
simulation box. As sketched in panel (a) of fig. [l] we 
perform plane shear transformations of the instantaneous 
shear strain 7 -A 7 + £7 with an essentially infinitesimal 
strain increment £7. We assume that not only the box 
shape is changed but that the particle positions r (us¬ 
ing the principal box convention m follow the imposed 
macroscopic constraint in an affine manner according to 


r x ^r x + J7 r y for |^-y| < 1 (46) 

with all other coordinates remaining unchanged. Albeit 
not strictly necessary for the demonstration of our key 
relations, we assume, moreover, that this shear transfor¬ 
mation is also canonical H EQ- This implies that the 
x-component of the velocity must transform as m 


V x v x — v y with |Jy| 1. 


(47) 


We emphasize the negative sign in eq. (47) which assures 
that Liouville’s theorem is obeyed OS ET 


C. Shear stress and affine shear-elasticity 


Let = Hid(£7) + Hex(£7) denote the system 

Hamiltonian of a configuration originally at 7 strained 
using a canonical affine transformation to 7 + £7 com¬ 
pactly written as a function of the strain increment £7. 
The instantaneous shear stress f and the instantaneous 
affine shear-elasticity Aa may be defined as the expansion 
coefficients associated to the energy change 

5H(5~i)/V = fSj + Aa^ 7 2 /2 for |Jy| 1 (48) 

with 7 being the reference, i.e. |32 

f = and (49) 

Aa m n"(5 7 )/v\^ (50) 


The derivatives H^Jy) and T-L"( 8 7) with respect to £7 
may be computed as shown in sect. 2.1 of m- Similar 
relations apply for the corresponding contributions rid 
and f ex to f = fid + f ex and for the contributions /byid 
and //A, ex to Aa = Aa, id + A a, ex- Using eq. (|44| this 


implies m 


Tid 


MA,id 


1 

V 


1 N 

i= 1 
N 


Yj miV lv 

i=l 


■Vi^y and 


(51) 

(52) 


for the ideal contributions to the shear stress and the 
affine shear-elasticity. Note that the minus sign for the 
shear stress is due to the minus sign in eq. (47) required 
for a canonical transformation. For the excess contribu¬ 
tions one obtains ITO] 


Tex = ^7 y7/i/(r-/) m }X ni } y and (53) 

A A.ex = y^{r?u"(ri)-rtu\n))nl x nl y 

+ ^Y^nu'in) n\ y (54) 


with = rjri being the normalized distance vector 
r = r_j — r r between the particles i and j. Interest¬ 
ingly, eq. ( [53] ) is strictly identical to the corresponding 
off-diagonal term of the Kirkwood stress tensor pQ. The 
last term in eq. (54) automatically takes into account the 
finite normal pressure of the system. 


D. Groundstate characterization 

Specific network. As explained elsewhere [8j [10] , the 
specific elastic network used in this work has been con¬ 
structed using the dynamical matrix of a quenched poly- 
disperse LJ bead glass, i.e. at low temperatures our net¬ 
work has exactly the same mechanical and vibrational 
properties as the original discrete particle system. Prior 
to forming the network the bead system was cooled down 
to T = 0 using a constant quenching rate and imposing 
a normal pressure P = 2. The original LJ beads are 
represented in the snapshot shown in fig. [2] by grey poly- 
disperse circles, the permanent spring network created 
from the quenched bead system by lines between vertices. 
The dark (black) lines indicate repulsive forces between 
the vertices, while the light (red) lines represent tensile 
forces. The line width is proportional to the tension (re¬ 
pulsion). Note that the force network is strongly inho¬ 
mogeneous with zones of weak attractive links embedded 
within a strong repulsive skeleton as already discussed 
in refs. [24l [25] . Only a small subvolume of the network 
is represented. The total periodic box of linear length 
L ~ 102.3 contains N = 10 4 vertices and Ni = 9956 
springs. The monomer density p is close to unity. 

Finite shear stress To . Due to the construction of the 
network the total force acting on each vertex of the ref¬ 
erence network must vanish at T = 0. As seen in the 
snapshot, this does not imply that the repulsive and/or 






0.20 



the shear stress r(7) as may be seen from the large cir¬ 
cles indicated in fig. [2j As indicated by the bold solid 
line these final stresses r(y) scale again linearly as 

t( 7 ) = to + G eq 7 = G eq (7 - 70) (55) 

with G eq ~ 16.3 and 70 = —To/G eq ~ —0.00072. The 
affine coefficient /xa has thus been replaced by the much 
smaller equilibrium shear modulus G eq of the ground- 
state network. Note that shear stress vanishes at a strain 
7 = 70 as indicated by the vertical arrow. If the strain 
7 is allowed to change freely (e.g., using a steepest de¬ 
scend scheme) without any external force r ex t applied, 
the system relaxes to 7 = 70. 


E. Finite temperature simulations 


FIG. 2 : Groundstate of elastic network model considered in 
this work assuming eq. ( |45| : (a) Snapshot of a small subvol¬ 
ume of linear length 10 containing about 100 vertices. The 
lines represent the quenched forces of the athermal (T = 0) 
reference configuration, (b) Shear stress r( 7 ) assuming an 
affine strain (squares) and after energy minimization (cir¬ 
cles). In agreement with ref. [10] we find in the first case 
t( 7 ) = to + /xa 7 with /xa = 34.2 (dashed line) and in the sec¬ 
ond case r(y) = To + G eq 7 with G eq = 16.3 (bold solid line). 
The stress does not vanish at 7 = 0 where r — To = 0.01161 
(horizontal arrow) but at 7 = 70 = —0.00071 (vertical arrow). 


tensile forces transmitted along the springs must also 
vanish. Due to the periodic boundary conditions and 
the constant-strain constraint (7 = 0) the shear stress r 
does not necessarily vanish. For a pair potential such as 
eq. (45) the relevant excess contribution f ex of the shear 
stress is readily computed using the Kirkwook expression, 
eq. ( [53] ) . As shown in the main panel of fig. [2] (horizon¬ 
tal arrow), it turns out that for the specific network we 
use throughout this work we have a finite shear stress 
To = r( 7 = 0) = 0.01161. 


Affine shear-elasticity /x a- Let us consider a small 
affine shear strain according to eq. ( [46] ) . As show in 
panel (b) of fig. [2] one may now compute using eq. (53) 
the shear stress r(q) for different 7 (squares). As one 
expects from the definition of the affine shear-elasticity 
coefficient /xa, eq. ( [48] ) , this yields the linear relation 
r(q) = To + p a 7 with /xa ~ 34.2 as indicated by the 
dashed line. The same coefficient /xa is also obtained di¬ 
rectly from the unstrained configuration using eq. (54). 


Equilibrium shear modulus G eq . The forces /. acting 
on the vertices i of an affinely strained network do not 
vanish in general and the system is normally not at me¬ 
chanical equilibrium. As described elsewhere [K^ [25], we 
relax these forces by first applying a steepest descend al¬ 
gorithm, i.e. by imposing displacements proportional to 
the force, and then by means of the conjugate gradient 
method [6j. The ensuing non-affine displacements of the 
vertices decrease the energy m and the magnitude of 


Molecular dynamics scheme. As discussed in sect. |IV[ 
this network is investigated numerically basically by 
means of a molecular dynamics (MD) simulation [lj [4] 
at constant particle number N = 10 4 , box volume 
V ~ 102.3 2 and a small, but finite mean temperature 
T = 0.001. Newton’s equations are integrated using a 
velocity-Verlet algorithm with a tiny time step ££md = 
10 -4 . The temperature T is fixed using a Langevin ther¬ 
mostat with a relatively large friction constant ( = 1. 
This was done to suppress artificial long-range correla¬ 
tions in the two-dimensional periodic simulation box. 

Thermal averages. Using the measured instantaneous 
shear stress f and affine shear-elasticity fi a we compute, 
e.g., the thermal averages t = (f), /x rr = /3V (Sr 2 ) and 
/xa = {ft a)- It can be shown for the respective ideal 
contributions that m 

Ed = Prr ,id = PA,id = F\d ~ ^ P (56) 

with Pid being the ideal normal pressure contribution. 
This holds irrespective of the considered A-ensemble. The 
ideal contributions are thus negligible. Due to the equiv¬ 
alence of the different axes it is also seen for the excess 
contributions that m 


MA, ex = Mb - Ax with 

MB = ~ r l U 'i r l)) n lx n ly) ( 57 ) 

l 


being the Born-Lame coefficient |7, I22H251 135] and P ex 
the excess part of the normal pressure P = Pid + P ex - 
Monte Carlo shear-barostat. The MD algorithm for 
the particles is coupled for A < 1 with a Monte Carlo 
(MC) scheme [T] [5] attempting every MD time step a 
canonical affine strain increment 7 -A 7 + ^7 (sect. |mB| . 
First, a strain increment £7 is randomly chosen from 
a uniformly distributed interval [-^7 ma x^7max]- In or¬ 
der to determine the Metropolis weight [5] we compute 


next the energy change STL = ?7(7 + £7) — TL{ 7) of 
the network which comprises both an excess contribu¬ 
tion due to eq. (46) and an ideal contribution due to 
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Gext 

A 

M77 

Ptt 

A 

77 

T 

j.77 

T 7r 

6W 

-10 

-1.59 

0.159 

60.1 

0.997 

0.0081 

9.2 

4.3 

13 

-6 

-0.58 

0.097 

43.7 

0.996 

0.0103 

5.6 

2.7 

6 

-3 

-0.23 

0.075 

37.9 

0.996 

0.0117 

4.3 

2.1 

4 

0 

0 

0.061 

34.2 

0.996 

0.0130 

3.5 

1.7 

2 

16.3 

0.50 

0.031 

26.1 

0.994 

0.0184 

1.8 

0.85 

1 

100 

0.86 

0.009 

20.1 

0.993 

0.0348 

0.5 

0.24 

0.3 

1000 

0.98 

0.001 

18.1 

0.979 

0.1030 

0.06 

0.03 

0.3 

10000 

0.998 

- 

17.9 

0.935 

0.3233 

0.007 

0.003 

0.4 


TABLE I: Summary of some properties as a function of the 
external spring constant G ex t- A = G ex t/(G eq + G ex t) with 
G e q = 16.3, strain-strain fluctuation /x 77 (fig. pi), stress-stress 
fluctuation p TT (fig. [7]), acceptance rate A of MC shear- 
barostat (fig. [3]), rj( A, /^-parameter defined by eq. (61), strain- 
strain crossover time T 77 ~ # 77 (fig# , strain-stress < 
time Ty r 

(fig- [is)- 


dry 


; crossover 

(fig. [lO]) , stress-stress relaxation time 0 TT 
The dynamical properties (columns 5-9) are only 
given iui ri — 10 -2 . The time scales T 77 , T 7T and 0 TT should 
be compared to the intrinsic relaxation tim e ta ~ 0.13 of the 
stress-stress correlations for n —>• 0 (fig. ED- 


eq. (47). Since our system is subjected to an external 
field E4xt(7), eq- (§, this gives rise to an additional con¬ 
tribution 5E4xt = C/ext (7 + ^7) ~ E/ext (7) • The suggested 
strain move £7 is accepted if 


£ < exp[—/3(<5"H + <5f/ ext )] 


(58) 


with £ being a uniformly distributed random variable 
with 0 < £ < 1 [5]. 

Operational parameters A and n. We assume for the 
external field E/ext (7) that r ext = 0 and 7 ext = 70, 
i.e. the average shear stress r is imposed to vanish 
for all ensembles studied and all ensembles compare the 
same thermodynamic state point. (This was explicitly 
checked.) The only remaining operational parameter 
from the static point of view is the external spring con¬ 
stant G ext or, equivalently, A = G ex t/(G e q + G ex t) as 
sketched in panel (c) of fig. [l] As seen in table [I] or 
fig- m we vary G ex t from —10, i.e. A ~ —1.59, over 
G ex t = A = 0 (NVrT-ensemble) up to G ex t = 10000, 
i.e. A « 0.998. The latter case is essentially equivalent 
to the standard NVyT-ensemble, i.e. the strain fluctua¬ 
tions become irrelevant for most properties. The second 
operational parameter of this study is the maximum at¬ 
tempted strain displacement ^7 max which determines the 
impact of the shear-barostat on the relaxation dynamics. 
Since the Metropolis MC move for the strain is performed 
every MD time step of length ^md, it is convenient to 
use instead of ^7max the maximum strain increment rate 
n m ^7max/^MD- (Since all simulations are performed 
with the same ££md> ^ 7max and n are strictly equivalent.) 
We compare below the dynamical strain and stress cor¬ 
relations for the five rates n = 1,10 _1 ,10 -2 ,10 -3 and 
10 -4 for a broad range of A. For A = 0.5 we have com¬ 
puted in addition the values n = 3, 0.3,0.03, 0.003 and 


1.0 


0.8 - 


0.6 


0.4- 


0.2 - 


0.0 


><0 | < 3 ) 48 X 3 ^ 

range of focus 


range c 
for dynamics 


X G ext =-10, >.=-1.59 
O G ex,=-9, A=-1-23 
□ G ext =-8, >.=-0.96 
> G ext =-6, >i=-0.58 
< G ext =-3^=-0.23 
OG ext =0A=0 
0 G ext =16.3, >1=0.50 
+ G ext =100, ?t=0.86 
A G ext =1000, >1=0.98 
V G ext =10000, A=0.998 




< 


+ 

good range 
for statics 


L. 


10 


10 


10 “ 


10 




10 


10 


10 “ 


FIG. 3: Acceptance rate A of the Metropolis MC shear- 
barostat, eq. (58), for a large range of A and n as a function 
of the scaling variable rj{ A, k) as defined by eq. (61). 


0.0003. Some properties for n = 10 -2 are summarized 
in the columns 5-9 of table [H The simulations become 
increasingly time consuming with decreasing k and data 
obtained for n = 10 -4 have to be taken with care. 

Free strain diffusion limit. The impact of k can be 
best judged from the simple limit where neither the sys¬ 
tem nor the external field restricts the barostat. This is 
the case for sufficiently small n depending on A. In this 
limit one expects the free diffusion of the strain y(t), i.e. 
g 77 (t) = Dryryt/2 for t <C T 77 . Since the (attempted and 
accepted) strain increment £7 is uniformly distributed 
in [—^7max, ^7max]? we have a mean-square strain step 
(^7 2 ) = ^7max/3 every ££md- This implies a diffusion 
constant 




pv 67 ; 


max 


3 StMB 




(59) 


in the free-diffusion limit. Using eq. (28) and eq. @ this 
yields in turn the corresponding crossover time 


T 77 = 6^y = 6 dtMD/rf 


y 2 

/ max 


pv 57* 

where we have defined the dimensionless variable 

n ^md 


rj(\, k) = yATTTT = 


AW) 


(60) 


(61) 


Since G eq and ^md are kept constant in all our simula¬ 
tions, the parameter 77 determines the dynamical regime 
for a system of operational parameters A and k. As seen 
in fig. [3] using 77 as a scaling variable the acceptance rate 
A of the MC shear-barostat of a broad range of A and 
n can be brought to collapse. For 1 <C 77 <C 10 the 
MC barostat is most efficient for static properties. For 
dynamical properties we shall focus below on /^-values 
where 77 <C 0.1 and thus A ^ 1 as emphasized by the 
solid horizontal line. 
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FIG. 4: Normalized histogram p{x) of x — 7 —7ext for different 
A. The line indicates for A = —1.59 (G ex t = — 10 ) a Gaussian 
with (Sj 2 ) s 1/pVGeft « 0 . 00012 2 . 



FIG. 6 : Strain-stress correlations: (a) Pre-averaged shear 
stress t(x) as a function of x — 7 — 7 ex t- All data collapse on 
the linear slope with G eq = 16.3 (bold solid line), (b) Scat¬ 
ter plot of 100 data points ( 7 , f) for A = —1.23 (G ex t = —9) 
with circles corresponding to a tiny time interval St = 10 -4 
(dashed line) and squares to a large St — 10 +2 (bold solid 
line) between subsequent configurations. 


The given data refer to a maximum strain rate k = 10“ 1 
for the Metropolis MC step 7 — )> 7 + £7 attem pted after 
each MD time step of %d = 10 -4 (sect. HIE). That this 
leads indeed to the desired Gaussian distributions can be 
seen from fig. [4] where the normalized histogram p(x) is 
plotted as a function of x = 7 — 7 ex t with 7 ext = 79. 
As one expects, p{x) S(x) for A —>► 1. The rescaled 
second moment /i 77 e /3V(S^/ 2 ) of the distribution is fur¬ 
ther analyzed in fig. |5j As one may see from panel (a), 
y( A) = /i 77 G eq decreases linearly with A (bold solid line) 
and vanishes for A = 1 as it should according to eq. ( [T()| ) . 
In turn this implies that one may determine G eq from 


FIG. 5: Second moment of strain fluctuations: (a) Dimen¬ 
sionless y “ /x 77 G eq with G eq = 16.3 as a function of A 
showing the linear decay (bold lin e) expected from eq. 0 - 
(b) Shear modulus G 77 , eq. ( [62] ) , confirming G eq ~ 16.3 
(bold line), (c) MSD g 11 {t) vs. time t for several A and 
k = IQ" 1 . T he dash-dotted line indicates the free-diffusion 
limit, eq. ( |59| ), for short times, the horizontal lines the long¬ 
time limit /x 77 for each A. Also indicated is for A = 0 and 
k = 1 (filled circles) an example with g 0.1 where the 
initial diffusive regime is suppressed. 


IV. COMPUTATIONAL RESULTS 
A. Static properties 

Strain-strain fluctuations. The equilibrium fluctua¬ 
tions of the strain 7 at a finite temperature T = 0.001 for 
different values of G ex t or, equivalently, A = G ex t/(G eq + 
G ex t) with G eq = 16.3 are presented in fig. [4] and fig. [5j 


G 77 — 1 /M 77 ^ext* (62) 

As can be seen from panel (b) of fig.|5j using this relation 
one confirms that the groundstate value G eq ~ 16.3 re¬ 
mains unchanged at a finite temperature T <C 1. The ac¬ 
curate determination of G 77 becomes of course more del¬ 
icate with increasing A since the MC acceptance rate be¬ 
comes eventually too small. The data quality in this limit 
can be readily improved (not shown) by either increasing 
the sampling time t tra j or by decreasing k, ~ rj. (The 
sampling becomes again inefficient if 77 gets too small.) 

Strain-stress correlations. Since the shear modulus 
G eq is finite, this implies that shear strain and shear 
stress fluctuations must be correlated 0. These corre¬ 
lations are addressed in fig. [6] The main panel presents 
the pre-averaged stress r(x) for an instantaneous strain 
x = 7 — 7 ext for several A [34]. Irrespective of x or A 
all data collapse on the linear slope indicated by the 
bold line with G eq = 16.3. Naturally, it follows from 
fig. [4] that the statistics must strongly decrease for x 2 
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X 


FIG. 7: Determination of shear modulus G eq ~ 16.3 (bold 
solid line) using G 77 for A < 1, G 7T for A < 1 and G TT for 
A ^ 0. Also indicated are the affine shear-elasticity /xa(A) ~ 
34.2 (dashed line) and the stress-stress fluctuation g TT which 
is seen to decay linearly (bold dash-dotted line) according to 
eq. ^ down to g T r\x=i ~ 18. 


fi 1T = /3V(S^Sr), i.e. a larger linear-slope window is vis¬ 
ible in fig. [6] for smaller A. Panel (b) of fig. [6] presents 
two scatter plots of (7,f) for A = —1.23 (G ex t = — 9). 
The small circles correspond to a time series of 100 subse¬ 
quent and, hence, strongly correlated configurations with 
St = StyiD . On these time scales the configuration has 
no time to relax the affine displacements imposed by the 
shear-barostat. As emphasized by the dashed slope this 
leads to a linear stress-strain relation with a coefficient 
/iA = 34.2. Other short-time sequences yield similar, 
but horizontally shifted linear slopes (not shown). A 
^-independent representation of the static strain-stress 
correlations for large times is obtained if (7, f) is indi¬ 
cated for 100 configurations with a time interval St = 10 2 
(squares). The scatter plot of the latter data is already in 
nice agreement with a coefficient G eq = 16.3 (bold line). 
Sampling over all data tuples (7, f) of a given A-ensemble 
one verifies that the shear modulus G eq is accurately de¬ 
termined using the linear regression coefficient 


V 7 / M77 


(63) 


Since, as shown in sect. |II A , (i 11 = (1 — A)/G eq and 


H 1T = 1 — A, this ratio must yield G eq for all A < 1. 
Confirming G eq = 16.3 the values of G 7r for different A 
are identical to G 77 as shown in fig. [7] 

Stress-stress fluctuations. The stress-stress fluctua¬ 
tions /i rr = /3V(St 2 ) presented in fig. [7] (circles) decrease 
linearly with A as stated in the Introduction, eq. ©• 
Being a simple mean /iA (small squares) is found to 
be strictly A-independent. As expected, fi TT -A /i a for 
A -A 0. Using eq. © the generalized stress-fluctuation 



FIG. 8: Diffusion of instantaneous strain 7 (t): (a) MSD 

g 11 {t) for A = 0 (open symbols) compared to the explicitly 
computed creep compliance J(t) = (S'yft)) / 5r ex t for a step 
stress 5Text — 0.01 (filled symbols). The dash-dotted lines in¬ 
dicate the regime used to determine the diffusion coefficient 
Z) 77 , the bold horizontal line the expected asymptotic limit 
/z 77 and the thin solid line the Kelvin-Voigt model, eq. (40), 
for n = 10 -2 . (b) Double-logarithmic representation of the 
rescaled crossover time T 77 /6 (Vmd vs. 77. The dash-dotted 
line indicates the expected power law, eq. (60). 


formula for G eq reads 


Gn-n- = 


A l^TT 


1 (/^A /^rr)/Gex t 


for A ^ 0. (64) 


For A —1, i.e. G ex t 00, this reduces to eq. (16). Equa¬ 


tion (64) thus generalizes the stress-fluctuation formula 


for the NVyT-ensemble [4j [8l[T0l I2TU231 [33] to general 
A. As seen from fig. [ 7 ] (diamonds), it is thus possible to 
determine G eq from fi a and fi TT for all A 7^ 0. 


B. Dynamics: Strain-strain correlations 

Strain-strain MSD. The strain histograms p(x) for 
different A and their second moments shown, respectively, 
in fig. [4] and fig. [5] have been rapidly sampled by aver¬ 
aging over short time series from trajectories of length 
t t raj = 10 4 . This is demonstrated in panel (c) of fig. [5] 
presenting the MSD g 17 (t) of the instantaneous strain 
7 (t) for k = 10 _1 . The statistics is improved by perform¬ 
ing a gliding average over the 10 8 data tuples sampled pQ. 
It is seen that p 77 (t) converges rapidly after a time T 77 
of order unity to its asymptotic limit /i 77 (A) indicated 
by horizontal lines. Free diffusion is observed for short 
times, i.e. g 7q {t) ~ D 77 t/2 for t <C T 77 (dash-dotted 
line), and this essentially irrespective of the ensemble in 
agreement with eq. ( [59] ) . Note that free strain diffusion 
is observed for all A and sufficiently small k, if the param¬ 
eter rj( A, k) < 1. As shown in the main panel of fig. [8] 
for different k and A = 0 (open symbols), this behavior 
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FIG. 9: Determination of equilibrium shear modulus G eq as a 
function of the sampling time t for G 11 {t) (thin lines), G 1T (t) 
(filled symbols) and G TT {t) (large open symbols) for k = 10 _1 . 
All fluctuation formulae are equivalent for t > 1 and converge 
similarly to G eq for large times. Note that G 7r (£) for all A < 1 
is everywhere identical to G TT (t) for A = 1. 


FIG. 10: Scaling of strain-stress MSD g 1T {t) for different A: 
(a) scaling collapse of y = g 1T {t)/ g 1T vs. x = t/T ir for n = 
10 -2 with the dash-dotted line indicating the diffusive regime 
and the bold solid line the long-time limit, (b) crossover time 
T 7r , eq. ([28]), rescaled as T 7T /65£md vs. the scaling variable 
r) using the same symbols as in the main panel. The dash- 
dotted line indicates the prediction eq. (65). 


may be used to determine first the short-time diffusion 
coefficient D 77 for systems with sufficiently small g and 
then using eq. ( [28] ) the crossover time T 77 of the strain 
fluctuations. The rescaled times T 77 /65 £md are repre¬ 
sented in panel (b) of fig. [8] as a function of the scaling 
variable g for several A. The dash-dotted line indicates 
the power-law slope, eq. ( |60| ). A perfect data collapse on 
this line is observed for all g <C 1. A successful data col¬ 
lapse can be also obtained for g 77 (t) over a broad range 
of A and k, by plotting # 77 (t)//i 77 as a function of the 
reduced time x = £/T 77 (not shown). As seen by the 
thin solid line in the main panel of fig. [8j one verifies 
that the scaling function f(x) = g 77 (t)//i 77 is given by 
f(x) = 1 — exp(— x). This is, of course, consistent with an 
exponentially decaying strain-strain relaxation function 
G 77 (£) = /i 77 exp(— x) as we have also checked directly. 

Creep compliance J(t). As discussed in sect. |fTe 1 
a MSD g 17 (t) computed in the NVrT-ensemble corre¬ 
sponds to a creep compliance J(t) = (Sj(t)) /Ar ext mea¬ 
sured by imposing at t = 0 a step stress increment 
15r e xt | <C 1 to the external shear stress r ext applied to the 
system. Using Ar ext = 0.01 and averaging over 100 config¬ 
urations this yields the data shown by the filled symbols 
in panel (a) of fig. [8] for two values of k. An excellent data 
collapse g 77 (t) « J(t) is observed. While it is natural to 
use the NVrT-ensemble for obtaining the creep compli¬ 
ance from the equilibrium fluctuations, it is worthwhile 
to note that due to the scaling ^ 77 (t)//x 77 = f(x), J(t ) 
may be also determined from other A-values. 


Time-dependent strain-strain fluctuations. As ex¬ 
pected from sect. |IIB| the strain-strain fluctuations 
J1 (t) computed according to eq. (22) by gliding aver¬ 
age from a finite time series increase monotonously with 


time t from zero (if only one configuration is sampled) to 
the asymptotic thermodynamic limit /i 77 (not shown). 
Note that one determines first the strain fluctuations in 
an interval [tofli = to + t\ and averages then over all 
possible to- If the estimate G 77 , eq. ([62|), of the modulus 


G eq is determined using /I (t) instead of /i 77 , it becomes 
a time dependent quantity called G 77 (t). As can be seen 
from fig. [9] for three different A values, G 77 (t) thus be¬ 
comes a monotonously decreasing function of time (thin 
lines). Interestingly, G 77 (t) appears not to depend on A. 
We note finally that using either the large-time asymp¬ 
totics eq. (27) or, more conveniently, the Debye relation 
eq. (29), one may determine the relaxation time # 77 of 
the strain-strain fluctuations J1 (t). One verifies that 
$ 77 T 77 holds as expected f or an exponentially decay¬ 

ing correlation function (sect. E5 


C. Dynamics: Strain-stress correlations 


Strain-stress MSD. The scaling of the strain-stress 

The relaxation time 


MSD g 1T (t) is investigated in fig. 10 


T 7r presented in panel (b) has been determined, as be¬ 
fore for T 77 , by matching the short-time diffusive regime 
with the long-time limit g 1T (t) —>> /i 7T . As indicated by 
the dash-dotted line, the data are well described by 


G ( 


= \j k 1 for 77 1 . 

Pa 


(65) 


That this scaling must hold can be seen by replacing 
in the definition of g 1T (t) the stress fluctuation Sf(t) by 
(t) which implies D lr = for the diffusion 
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coefficients and in turn eq. (65) using again eq. ( [28] ) . We 
stress that in agreement with the inset of fig. [6] the system 
has not enough time to relax the affine strain imposed by 
the shear-barostat for the short times t <C ta ~ 0.1 con¬ 
sidered here. It is thus g A and not G eq which has to be 
used as the linear slope coefficient. The main panel (a) 
presents the collapse of y = g 1T (t)/g 1T as a function of 
x = t/T lr for ft = 10 -2 . Minor deviations are visible for 
x « 1. The quality of the collapse improves by decreas¬ 
ing 77(A, ft), whereas the crossover becomes more sudden 
(even a hump may occur) if 77 is too large (not shown). 

Time-dependent correlation coefficient G 1T {t). Since 
the strain-shear correlations are dominated for t < r A 
by the affine shearing, one expects G lr (t) computed us¬ 
ing the ratio of ~jJ 1T (t) and ~p, (t) to be similar to g A 
for t —>• 0. As one sees from G lr (t) presented in fig. [9] 
(filled symbols), this is indeed the case. Interestingly, 
G yr (i) is seen to be independent of A for all times £, i.e. 
the A-dependences of ~p 7T {t) and ~pL (t) do cancel for all 
times just as they cancel for the ratio of the equilibrium 
moments /i 7r //x 77 . 


D. Dynamics: Stress-stress correlations 

Introduction. We turn finally to the more intricate 
characterization of the stress-stress correlations as a func¬ 
tion of our two operational parameters A and ft. We focus 
again on the regime with 77(A, ft) 1. As above we begin 
by discussing the MSD g TT (t). Since according to eq. (|20|) 


g TT (t) = C TT ( 0) - C TT (t ) with C TT { 0) = fi TT , (66) 

the MSD g TT (t ) and the correlation function C TT (t ) con- 
tain in principal the same information. From the presen¬ 
tational point of view g TT (t) has the advantage that the 
short-time power-law behavior of the correlations can be 
made manifest using a double-logarithmic plot as shown 
in fig. [Tl] We describe then the large time behavior of 
the correlation function C TT (t) and demonstrate that r* 
is given by T 77 . Finally, we turn to the time-dependent 
stress fluctuation ~p TT (t) which is used to determine the 
stress-stress relaxation time 0 TT . 

Stress-stress MSD. Let us concentrate first on the 
data for A = —1.59 presented by open symbols in fig. [TT] 
If ft is large, all internal dynamics is destroyed and we 
find immediately g TT (t) ~ g TT as seen for ft = 1 (small 
circles). For smaller, but not too small ft one observes 
a free-diffusion regime with g TT (t) = D rr t/2 as shown 
by the bold d ash-do tted lines on the left. As above for 
g 1T (t) in sect. IV C| Sr(t) in the definition of g TT (t) can 
be replaced by pL A o^{t). This argument yields the dif¬ 
fusion coefficient D rr = successfully used in the 

plot for ft = 10“ 1 and ft = 10 -2 . If ft is further reduced, 
g TT {t) becomes ft-independent for short times. As seen 
for ft = 10 -3 and 10 -4 the dynamics becomes “determin¬ 
istic” in 



FIG. 11: Stress-stress MSD g TT (t) focusing on A = —1.59 
comparing different ft (open symbols). Additionally, we 
present A = 0 for ft = 1 (small filled circles) and ft = 10 -4 , 
A = 0.5 for ft = 10 -4 and the NVyT-ensemble (stars). The 
dot-dashed lines indicate the free diffusive behavior for inter¬ 
mediate 77 , the thin solid line the analytic short-time limit for 
small 77 , the horizontal dashed line the thermodynamic limit 
g TT for A = —1.59 and the bold solid line the expected inter¬ 
mediate plateau, eq. (38). Note that for sufficiently small ft 


all data approach the ft- and A-independent lower envelope in¬ 
dicated by the solid li nes in agreement with the fundamental 
scaling postulate eq. (31). 


with 6^ ~ 0.16 as shown by the thin solid line. Albeit fi TT 
depends on A, the ratio /a TT /0 2 does not, in agreement 
with the fundamental postulate eq. (31). This implies 
that 6 depends somewhat on A. The t 2 -scaling in eq. (67) 
is expected [10] since the associated correlation function 
C TT (t) must be a continuous and symmetric function at 
t — 0 which, moreover, should be an analytic function if 
barostat and thermostat effects can be ignored [35]. This 
implies that g TT (t) must be an even expansion in terms 
of t 2 which to leading order leads to eq. (67). Since there 


9r.it) = (t/6Y 


~ ft u A u for t <C ta (67) 


are several dynamical regimes it is not meaningful to de¬ 
termine a crossover time T rr using eq. ( |28| ) as before for 
T 77 and T 7r over the full range of A and ft. We shall see 
at the end of this section how a longest relaxation time 
6 rr for the stress-stress correlations might be obtained. 
For times t > ta an intermediate plateau appears (bold 
solid line) as predicted by eq. (37). As can be seen for 
four different A values this intermediate plateau does not 
depend on the ensemble provided that ft is sufficiently 
weak. (As one expects from fig. [3] deviations from this 
asymptotic limit arise again for large ft if 77(A, ft) 0.1.) 
We emphasize that the two solid lines indicated in fig. [TT] 
give the lower envelope for all parameters A and ft we 
have simulated. The MSD g TT (t) behaves thus as an 
ensemble-independent simple average as we have argued 
in sect. m One may operationally define the crossover 
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FIG. 12 : Relaxation modulus G(t) obtained from the stress 
response to an applied strain £7 = 0.01 (stars) compared to 
the equilibrium correlation function C TT (t) | A for several A ob¬ 
tained by averaging over 1000 quenched-strain configurations. 



FIG. 13 : Effect of the barostat parameter ft for C TT {t) with 
A = 0 . 5 . The large stars correspond to the quenched 7- 
ensemble (r* = 00) and all other symbols to a switched- 
on barostat of finite ft. The vertical arrows indicate T 77 
for ft — 10 -1 ,10 -2 ,10 -3 and 10 -4 confirming that roughly 
T 77 ~ r*. A much longer run is warranted for k = 10 4 . 


time ta by matching eq. (67) and eq. (38). This yields 


r A = 6 


1 G e q/ MA 

1 AG eq / /Jj A 


0.13 


( 68 ) 


as indicated by the vertical line. Note that r A = 0V 2 
for A = 1. We stress that while 0 is a function of A, 
ta is an ensemble-independent constant for the given 
elastic network. (Moreover, it can be shown that it 
barely depends on the friction constant £ of the Langevin 
thermostat [ID].) Using eq. (J68| one may reformulate 
eq. ( [67] ) in a manifest A-independent form: g TT (t) = 
(/jl a — G eq ) (t / r A ) 2 . The MSD ultimately approaches 
/i rr (A) for even larger times t r* as indicated by the 
dashed horizontal line for A = —1.59. Remembering that 
fi TT =- fi A — G eq for A = 1 it is clear that the NVyT-data 
stays constant for all times. 

Relaxation modulus. The large-t scaling is further ad¬ 
dressed in fig. [l2] and fig. [13] presenting the stress-stress 
correlation function C rr (t). The most central result of 
this work stated by eq. m is demonstrated in fig. [l2j 
We compare the directly measured relaxation modulus 
G(t) with the equilibrium correlation function C TT (t) as¬ 
suming an asymptotically slow barostat (r* = 00). Av¬ 
eraging over 1000 independent configurations from an 
NVyT-ensemble with 7 = 0 for t < 0, the relaxation 
modulus G(t) = (f(t))/Sj measures the stress response 
after a (canonical and affine) shear strain £7 = 0.01 was 
been applied at t = 0 (stars). The correlation functions 
C TT (t) for A < 1 have been obtained by averaging over 
1000 equilibrated configurations from a A-ensemble as in¬ 
dicated. Switching off the shear-barostat the strain 7 of 
each configuration is quenched. Note that C TT (0) = fi TT 
holds due to the ensemble averaging. As expected from 
eq. (37), C TT (t) decreases monotonously from fi TT down 


to (1 — A)G eq for t r A . This corresponds to the 
A-independent intermediate plateau g TT (t) = fi A — G eq 
in fig. [TTJ It is seen that G(t) = C TT (t) for A = 0 
(large circles) and that for all A one may obtain G(t) by 
vertically shifting the correlation functions C TT (t )| A -A 
C rr (t) | A + AG eq = G(t) using the already determined 
shear modulus G eq . Note also that C TT (t) |a=i —» 0 for 
large times while all other C TT (t) remain finite. As al¬ 
ready emphasized at the end of sect. EH it is thus impos¬ 
sible to obtain the modulus G eq solely from C rr (t) |a=i- 
Please note that the observed short-time value fi TT « 18 
(bold dashed line) for A = 1 is quite different from the 
equilibrium modulus G eq ~ 16 (bold solid line). 

Finite-n effects. Focusing on one ensemble with A = 
0.5 we investigate in fig. [13] the effect of a switched-on 
shear barostat. The stars refer to data already seen 


in fig. 12 using the quenched-strain ensemble (k = 0, 
r* = 00). As one expects, the decay of C TT (t) be¬ 
comes systematically slower with decreasing k. For large 
k > 10 -2 we have computed C TT {t) by gliding averaging 
over a single trajectory of £ tra j = 10 4 - For smaller k it 
becomes increasingly harder to have a sufficiently long 
trajectory of {77J probing the phase space, i.e. for which 
G rr (0) = /i rr , and to store at the same time the shear 
stresses with a sufficiently fine time resolution. The data 
for /^ < 10“ 2 have thus been obtained by using an ensem¬ 
ble of 100 configurations of different 7 and averaging over 
the trajectories obtained for each configuration using a 
finite k. As expected from eq. ( [37] ) , for sufficiently small 
a all correlation functions have an intermediate plateau 
for ta < t C a (bold solid line). That the barostat is 
not completely switched off can be seen from the final de¬ 
cay. The vertical arrows indicate the crossover times T 77 
for the four smallest ft. As one expects, this time scale 
appears to coincide roughly (up to a constant prefactor) 
with the upper limit r* of the respective plateau. In other 
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FIG. 14: Shear-stress fluctuation ~p TT (t) | A as a function of the 
sampling time. Also indicated is the affine shear-elasticity 
fiA obtained for the NVrT-ensemble (small filled squares) 
which is seen to immediately converge to the long-time limit. 
The thin lines correspond to the integrated she ar-s tress auto¬ 
correlation function C TT {t) | A according to eq. (25). 


FIG. 15: Stress-stress relaxation time 0 TT : (a) Determina¬ 
tion of 0 TT from the decay of y — 1 — ~p TT (t) //jl tt accord¬ 
ing to eq. (27) for A = 0.5 for different ft. The thin dashed 


words, the plateau disappears when the ergodicity break¬ 
ing associated with the averaging over a quenched-strain 
ensemble is lifted by the shear-barostat. 

Time-dependent stress-fluctuation formula. Let us 
return to fig. [9| According to the generalized stress- 
fluctuation formula G rr , eq. (64), the shear modulus G eq 
may be obtained from the afnne shear-elasticity /ia and 
the equilibrium stress-stress fluctuation fi TT for A < 1 
and A ^ 0 (fig. [ 7 ]). If instead of y TT a time-averaged 
fluctuation ~p TT {t) is used, this leads again to a time- 
dependent monotonously decreasing estimate G rr (t) as 
seen from the open symbols in fig. [9j No shear-barostat 
is of course applied for the NVyT-ensemble indicated by 
the large squares. We emphasize the remarkable finding 
that within numerical accuracy 


lines indicate eq. ([27]), the bold horizontal line the expected 
intermediate plateau (1 — A)G eq / y T r for 10ta t <C t*. 
(b) Relaxation time 0 TT for three not too large A plotted as 
0 TT /65tMD vs. 7 y(A, ft). The dash-dotted line ind icates the 
power law expected for 77 <C 1 according to eq. (71). 


have integrated C TT (t) as suggested by eq. ( [25] ) . As seen 
by the thin lines indicated in fig. [l4j this leads to con¬ 
sistent results. If one has thus characterized G rr (£), this 
implies ~fl rr {t) and then in turn G rr (t). Note that for all 
A the asymptotic limit is reached at about IOOta which 


is of the same order as T, 


77 


3.5 for A = 0 and ft = 10 


G tt (£)|a=i ~ G qr (t) |a< 1 


(69) 


for all times and a broad range of ft for G ir (t). As seen for 
A = —0.58, for small A and not too short times G rr (t) and 
G 77 (£) are generally similar. Since ~p TT (t) must vanish for discussion in sect. 


small times, this implies that G TT (t) must become finite 


Relaxation time 0 TT . The convergence of Jt rr ( t ) is sys¬ 
tematically analyzed in fig. [15] The main panel presents 
the dimensionless deviation y = 1 —~p TT (t) / fi TT for several 
ft and A = 0.5. As emphasized by the thin dashed lines, 
a correlation time 0 TT may be measured from the large¬ 
time 1/t-decay following eq. (27). (Unfortunately, the 
sampling time is yet insufficient for our smallest ft.) With 
decreasing ft the decay becomes systematically slower 
and a horizontal plateau with y = (1 — A)G eq / pb TT (bold 
solid line) appears in an intermediate time window. This 
plateau is expected from eq. (30) and eq. ( |37| ). Following 
II B[ this implies that for suffi- 


G TT (fl) 


RA 


1 M A / G e xt 


for t -A 0 


(70) 


ciently small ft and not too large A the relaxation time 
0 TT must scale (up to a prefactor of order unity) as 


and, hence, G rr (t) -A fi a for G ex t Ma- This limit is 
indicated by the bold dashed line for A = 1. 

Time-dependence of stress-stress fluctuations. The 
time-dependence of ~fi TT (t) for several A is explicitly rep¬ 
resented in fig. 14 assuming ft = 10“ 1 for A < 1. While 


(1 - A)G ( 


eq 


the simple mean /ia converges immediately with time, 
~p TT (t) increases monotonously from zero to the large¬ 
time limit fi TT (horizontal lines). As already pointed out 
in sect. EH assuming time-translational invariance the 
time-dependence of ~p TT (t) is a consequence of the time- 
dependence of the C rr (t). To emphasize this point we 


r* ~ (1 — A) 2 /ft 2 for 77 1 (71) 

/XtT 

to leading order where we have used eq. ( |6Q| ) for r* ~ 
T 77 ~ 1/ft 2 . That the ft-scaling holds can be seen in 
panel (b) of fig. 

0 , 
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(For A = 0.5 we have computed 
also for the additional values ft = 3, 0.3, 0.03, 0.003 
and 0.0003.) The same representation as in the insets 
of figs. [8] and [lO] has been chosen for comparison. It 
should be stressed, however, that for 0 TT this is for two 
reasons not a rigorous scaling plot over the full range 
of operational parameters considered. Firstly, due to the 
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(1 — A )G eq /fi TT -factor in eq. (71) there is an additional A- 
dependence even for g <C 1 and not too large A. (This can 
be accounted for using a different representation.) More 
importantly, for A ~ 1 barostat effects naturally become 
irrelevant and the intrinsic relaxation time ta essentially 
sets 6 tt and not r* ~ T 77 . We thus observe 0 TT ~ taA 0 ^ 0 
for all A > 0.8 as seen in table [i] for n = 10 1 2 3 4 . Basically, 
since we have two time scales ta and T 77 , it is not possible 
to put all 0 TT ( A, k) on one master curve. 


V. CONCLUSION 


Static properties. Focusing on two-dimensional elas¬ 
tic networks (sect. Ill) we have characterized theoreti¬ 
cally (sect. [n| and numerically (sect. |IV[ ) the static and 
dynamical fluctuations of shear-strain and shear-stress 
in generalized Gaussian ensembles as a function of the 
dimensionless parameter A = G e xt/(G e q + G ex t), fig- l] 
Monitoring various properties we interpolate between the 
NVyT-ensemble (A = 1) and the NVrT-ensemble (A = 0) 
and consider even negative A. The static strain-strain 
fluctuations /i 77 (fig. [5), the strain-stress fluctuations g 1T 
and the stress-stress fluctuations g rr (fig. [7]) have been 
shown to decrease linearly with increasing A [20]. As a 
consequence, the static shear modulus G eq may be ob¬ 
tained (fig. [7]) from either the strain-strain fluctuation 
formula G^ 7 , eq. (62), or from the strain-stress formula 


G 7r , eq. ([63]) , for A < 1 or the stress-stress fluctuation 


formula G77, eq. (64), for A 7^ 0. The latter formula 


is a generalization of the well-known stress-fluctuation 
formula eq. ( [l6| ) for A = 1. When sampled from a fi¬ 
nite time series (7^,4), G 77 (t), G 7T (t) and G rr (t) be¬ 
have similarly with time t for all A (fig. 9|. Interestingly, 
G TT (t)\ x =i and G 7 T (£)|a<i are identical for all times. 

Dynamic properties. The infl uence of the parameter 
n of the MC shear-barostat (sect. |FQE ) has been investi¬ 
gated for various dynamical properties, especially for the 
stress-stress correlation function C TT (t) (fig.[l3|). The up¬ 
per time limit r*(A), below which barostat effects can be 
ignored, is set by the strain-strain crossover time T 77 (A) 
(fig. [8|. While a rapid barostat may be of advantage 
for static properties, a sufficiently slow barostat with 
a large r* corresponds to the fundamentally important 
linear-response limit where the barostat is only relevant 
for the distribution of the start points of the trajectories 


(and, hence, G a b(0) = /i a b) but not for their dynamical 
pathways for t < r*. We have argued that this implies 
the fundamental scaling g a b(£) ^ A 0 ft 0 for t <C r*(A) as 
demonstrated explicitly for g TT (t) in fig. [TT] Assuming 
such a slow barostat one may obtain Geq using eq. ([6| 
from C TT (t)\ x « (1 — A)G eq for A < 1 (fig.[l2|. As empha¬ 
sized elsewhere mm, it is impossible, however, to obtain 
the modulus G eq alone from the autocorrelation function 
G rT (£)| A =i- The experimentally important shear-stress 
relaxation modulus G(t ) may be most readily determined 
at A = 1 calculating first G eq = g a ~ hrr |a=i and then 
G{t) = G rr (£)| a=iTG eq [36]. We have commented briefly 
on the creep compliance J{t). Being well-described by 
the Kelvin-Voigt model, eq. ( |4Q| ), J(t) is best obtained 
from the strain-strain MSD g^fit) at A = 0 (fig. [8]). 

More general thermodynamic variables. Most of the 
theoretical relations and numerical techniques discussed 
above, especially the key formulae eq. © and eq. ©, 
generalize readily for any pair of continuous extensive and 
intensive variables X and I with M eq = dl/dX being the 
equilibrium modulus. With M ext being the spring con¬ 
stant of the external potential controlling the fluctuations 
of the extensive variable and A = M ex t/(M eq + M ex t) one 
obtains, e.g., the relaxation modulus 

M{t)=C(t)\ x=0 = G(t)| A + AM eq for t « t*(A) (72) 

with C(t ) = /3V(Si(t)SI( 0)) being the relevant au¬ 
tocorrelation function. The associated MSD g(t) = 
(/3V/2)((I(t) — /(0)) 2 ) is expected to be strictly A- 
independent in the same time window. These properties 
must be computed again either using a slow “barostat” 
with a large, albeit finite r* <C t tra j or, equivalently, by 
averaging over an equilibrium ensemble of configurations 
with frozen X. For times ^ > r* a switched-on “baro¬ 
stat” ultimately restores the ergodicity and C(t ) must 
vanish while g{t) approaches its finite thermodynamic 
value fiV(SI 2 ). 
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